CONVERGENCE OF A FINITE VOLUME SCHEME FOR GAS 
WATER FLOW IN A MULTI-DIMENSIONAL POROUS MEDIA 



MOSTAFA BENDAHMANE^ ZIAD KHALIL^, MAZEN SAAD^ 

^ Institut Mathematiques de Bordeaux, Universite Victor Segalen Bordeaux 2, Place de la Victoire, 33076 

Bordeaux, France. 

^ Ecole Centrale de Nantes, Universite de Nantes, Laboratoire de Mathematiques Jean Leray, UMR CNRS 
O ' 6629, 1, rue de la Noe, 44321 Nantes France 

On) , E-mail :mostafa. bendahmane@u-bordeaux2.fr, Ziad.Khalil@ec-nantes.fr, Mazen.Saad@ec-nantes.fr 

D . Research partially supported by GNR MOMAS 

(N 

^ ■ Abstract. A classical model for water-gas flows in porous media is considered. The degenerate 
1^ ■ coupled system of equations obtained by mass conservation is usually approximated by finite 
. volume schemes in the oil reservoir simulations. The convergence properties of these schemes 
J3 ! are only known for incompressible fluids. This chapter deals with construction and convergence 
analysis of a finite volume scheme for compressible and immiscible flow in porous media. In 
comparison with incompressible fluid, compressible fluids requires more powerful techniques. 
^ We present a new result of convergence in a two or three dimensional porous medium and under 

^ . the only modification that the density of gas depends on global pressure. 
OO 

o 



■ 1 Introduction 

o 



X 



Mathematical study of a petroleum-engineering schemes takes an important place in oil recov- 
ery engineering for production of hydrocarbons from petroleum reservoirs. In soil mechanics, 
engineers study the air-water flow in soils and they prefer the use of a two phase flow model. 



^ , More recently, due to the effects of global warming on climate change, two and multi-phase 
flow has been receiving an increasing attention in connection with the disposal of radioactive 
waste and sequestration of CO2. 

It has been shown that the governing equations describing two incompressible (compress- 
ible) phase flow in porous media can be written in a fractional flow formulation, i.e., in terms 
of global pressure and saturation and that formulation has been studied; For incompressible 
flow and from a mathematical point of view [?, [23] and it has been used in numerical codes 
[2T| [22l [20] . For immiscible and compressible two-phase flows (e.g., air, water), Ewing and 
al. in ([28], [25]) follow the ideas of Chavent by considering global pressure and saturation 
as unknowns of the system. This formulation leads to a global pressure equation coupled to 
the water saturation equation. The authors proposed a finite element and finite difference 
method to solve the saturation equation and a mixed finite element to approximate the global 
pressure equation. Note that the global pressure reads as a parabolic equation with a source 
term involving the evolution of the capillary pressure term. This evolution term is approached 
by Picard iterations. This algorithm converges numerically and suggests a continuous depen- 
dence on the capillary terms and legitimates some approximations for small capillary pressure. 
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Further, it has been proven that this fractional flow approach is far more efficient than the 
original two-pressure approach from the computational point of view [21] and the references 
cited therein. For compressible flow and from mathematical point of view, the fractional flow 
formulation is sufficient enough at least for slightly compressible gas, i.e, when the density of 
gas depends on the global pressure |32l [33] . More recently and under the context of theoretical 
study of compressible flow in porous media, the two-pressure approach has been treated by Z. 
Khalil, M. Saad [351 [36]. 

In this paper, we consider immiscible two-phase flows; the gas phase is considered to be com- 
pressible and the water one to be incompressible. The model is derived by using the global 
pressure notation and is justified at least for slightly compressible gas. The system represents 
two kinds of degeneracy. The first one is the classical degeneracy of the diffusion operator in 
saturation due to the capillary effect. The second one represents a degeneracy in the evolu- 
tion term in pressure occurring in the region where the gas saturation vanishes: A classical 
compactness result on pressure is missed in the region where the gas phase is missing. 

The aim of the present paper is to show that the approximate solutions obtained with the 
proposed upwind finite volume scheme f l3.5l) -( l3r6l) . converges as the mesh size tends to zero, to 
a solution of system (I2.ip - fl2.2p in an appropriate sense defined in Section |2l In Section [3] we 
introduce the finite volume discretization, the numerical scheme and state the main convergence 
results. In Section HJ maximum principle on saturation is attained and a priori estimates on 
the discrete gradient of the capillary term and on the discrete gradient of the global pressure are 
derived as the continuous case in C. Galusinski, M. Saad [33]- In Section a well posed- ness 
of the scheme is inspired by H. W. Alt, S. Luckhaus[?]. Section [6] is devoted to a space-time 

compactness argument, in this section we follow B. Andreianov, M. Bendahmane and R. 
Ruiz-Baier []. Finally, the passage to the limit on the scheme needs a powerful techniques due 
to the lack of compactness result on global pressure in the region where the saturation of gas 
vanishes, and this performed in section [71 

2 The mathematical formulation 

The fractional flow formulation describing the immiscible displacement of two compressible and 
incompressible fluids are given by the following mass conservation of each phase [33] : 




dti^s) + div(KM2(s)Vp) - div(Ka(s)Vs) + div(Kp2M2(s)g) + sfp = fp-fj. (2.2) 



where and K are the porosity and absolute permeability of the porous medium; p,p2,P 
and s are respectively the densities of gas and water (density of water is constant), the global 
pressure and the saturation of gas; fp,fj,Mi,M2 and g are respectively the production and 
injection source terms, the mobilities of gas and water and the gravity term. 
To define the capillary term a, let us denote by pi,P2 to be respectively the pressures of gas 
and water phases. Thus, we define the capillary pressure and the total mobility as 




dt{(t>p{v)s) - div(Kp(p)Mi(s)Vp) - div(Kp(p)a(s)Vs) 






(2.3) 
(2.4) 
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and the function s h-)- puis) is non- decreasing (^^(s) > 0, for all s G [0, 1]). In this paper, the 
forced displacement of fluids is modellized. It is used in many enhanced recovery processes: a 
fluid, such as water, is injected into some wells in a reservoir while the resident hydrocarbons 
are produced from other wells. Now we define the capillary term 

Mis) -dT^'^-' 

defining a function p{s) such that ^(s) = and setting p = P2+P, named global 

pressure [23]. Thus, each phase velocity given by Darcy's law can be written as 

Vi = -KMi(s)Vp - Ka(s)Vs + KMi(s)pi(p)g (2.5) 
V2 = -KM2(s)Vp + Ka{s)Vs + KM2{s)p2g. (2.6) 

Note that this system is strongly degenerate. In fact, the lack of coercivity of the degenerate 
diffusion term div(Ka(s) Vs) is classical for incompressible flows. An additional difficulty is due 
to the degeneracy of the time derivative term (j){x)dt{p{p)s) which vanishes in the region where 
s = 0. Another difficulty seems to be the degenerate diffusive pressure term div(Kp(p)Mi(s)Vp) 
in (12. ip where the mobility of the gas phase Mi vanishes in s = 0. In fact, a pressure diffusion 
term appears also in the saturation equation (12. 2 p with the term div(KM2(s)Vp). An energy 
estimate coupling the two equations (I2.ip - (l2.2p lets appear a non-degenerate dissipative pressure 
term (see section H]). 

Consider a fixed time T > and let Q he a. bounded set of M.'^ {d > 1). We set Qt = 
(0,T) X fi, St = (0,T) x dn. To the system ([53])-(E2]) we add the following mixed boundary 
conditions and initial conditions. We consider the boundary dQ = U Fj, where denotes 
the water injection boundary and Fj the impervious one. 



s{t,x) = 0, p(t,x) = on F^ 
Vi ■ n = V2 ■ n = on Fj, 

where n is the outward normal to the boundary Fj. We force a constant pressure (shifted at 
zero) along the time on the region of water injection. 
Initial condition: 

s(0, x) = Sn(x), in Q 
p(0, x) = po{x) in n 

We are going to construct a finite volume scheme on orthogonal admissible mesh, we treat 
here the case where 

K = kld 

where A; is a constant positive. For clarity, we take k = 1 which equivalent to change the scale 
in time. In remark mi we give the scheme where k is a function depending on space. 

Next we introduce some physically relevant assumptions on the coefficients of the system. 



(HI) The porosity (j) belongs to L°°{Q) and there exist two positive constants (/)o and such 
that 00 < (pix) < (pi a.e. x ^ Q. 
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(H2) The functions Mi and Ms belong to C°([0, 1]; M+), Mi(0) = and M2(l) = 0. In addition, 
there is a positive constant mo such that, for all s G [0, 1], 

Mi(s) + M2(s) > mo. 

(H3) The function a G ^^([0, 1]; M+) satisfies > for < s < 1, and a(0) = 0. 

We define (3{s) = a{z)dz and assume that /3~^ is an Holder function of order 9, with 
< 6' < 1 on [0,/3(l)]. This means that there exists a positive c such that for all 
Si,S2 e [0,/3(l)], one has \p-\si) - p-\si)\ < c|si - 

(H4) ifpjj) e {L\Qt))\ fp{t,x), fj{t,x) > a.e. {t,x) e Qt 

(H5) The density p is a C^(M) function, increasing, and there exist pm > 0, pm < +oo such 
that 

Pm < p{p) < Pm for all p eR. 

The assumptions (HlT])-(Hl5]) are classical for porous media. Especially, a practical sufficient 
condition to handle (HS]) is to consider that a is an Holder function at s = 0. This contains 
several relevant physical cases of two-phase flows in porous media (see [231 chapter V]). 

Define 

H^jn) = {ue H\n); n = on r^}, 

this is an Hilbert space when equipped with the norm (q) = ll^'"ll(L2{n))'*- 

In the next section we introduce the existence of solutions to system fl2.ip -( l2l2|) under the 
conditions (HUD -(HE]). 

Definition 2.1. (Weak solutions) 

Let (I{J^-(I{5^ hold. Assume Pq (defined by i\2.8^ ) belongs to L'^{Q), and Sq satisfies < So < 1 
almost everywhere in Q. Then, the pair {s,p) is a weak solution of the problem ( f^. jj) -( f^!^) if 

0<s<l a.e. inQT, (3{s) e L\0,T; H'^JQ)), p e L\0,T- H^JQ)), 

such that for allip,^ eV{[0,T)xn), 

(f)p{p)sdtip dxdt — / (f){x)uo{x)ip{0, x) dx 



+ I p{p)Mi{s)Vp-V(pdxdt+ / p{p)Vl3{s) -Vifdxdt 

Jqt 

p^{p)Mi{s)g-V'{)dxdt+ I p{p)sfp'^dxdt = 0, (2.9) 



n 



(psdtidxdt— I (f)So{x)S,(S^,x) dx + / V l3{s) ■ dxdt 

Jqt 

M2{s)Vp-Vidxdt- I p2M2{s)g-V^dxdt 



+ [ sfpidxdt= ! {fp- f,)idxdt. (2.10) 

JQt JQt 
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3 The finite volume scheme 

Now, we want to discretize the problem (12. ip - (12. 2p . Let T be a regular and admissible mesh of 
the domain Q, constituting of open and convex polygons called control volumes with maximum 
size (diameter) h. 

We let fl be an open bounded polygonal connected subset of with boundary dfl. Let T 
be an admissible mesh of the domain Q consisting of open and convex polygons called control 
volumes with maximum size (diameter) h. For all K & T, let by xk denote the center of K, 
N{K) the set of the neighbors of K i.e. the set of cells of T which have a common interface 
with K, by Ni^tiK) the set of the neighbors of K located in the interior of T, by Ne:>^t{K) the 
set of edges of K on the boundary dfl. 

Furthermore, for all L G Nint{K) denote by dK,L the distance between xk and xl, by (Tk,l the 
interface between K and L, by r]K,L the unit normal vector to aK,L outward to K. And for all 
a G Ncxt{K), denoted by dK,a the distance from xk to a. 




Figure 1: Control volumes, centers and diamonds 

For all G T, we denote by \K\ the measure of K. The admissibility of T implies that 
Q = U/^gT-i^, K n L = ^ \i K,L G T and K ^ L, and there exist a finite sequence of 
points {xK)K^r and the straight line xkXl is orthogonal to the edge aK,L- We also need some 
regularity on the mesh: 

mm — — — > a 

K&T,L(^N{K) diam(A ) 

for some a G IR'*'. 
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We denote by Hh{Vi) C L'^{Vl) the space of functions which are piecewise constant on each 
control volume K e T. For all Ufi e Hfi{Q,) and for all K e T, we denote by Uk the constant 
value of Ufi in K. For {uh,Vh) e (if^(Q))^, we define the following inner product: 



KeTLeN{K) ^'^ 

In the case of homogeneous Neumann boundary condition, for example Vu ■ rj = Vv • = on 
Fj C dfl, so we impose Ui^ — Uk = — Vk = H (Tk,l C Fj. And in the case of homogeneous 
Dirichlet boundary condition u — v — Q onV^j C dD,, so we impose Ul — Vl — ii (7k,l C F^„ 
and dK,L denotes the distance form xk to (Tk,l-i more precisely, 

We define a norm in Hfi{il) by 

Finally, we define Lh{0,) C L'^{Q) the space of functions which are piecewise constant on each 
control volume K & T with the associated norm 

for {uh,vh) e {Lhin)y. 

Next, we let X e T and L e N{K) with common vertexes (a^,ft:,z,)i<^</ with I e n*. Next 
let Tk^l (respectivley T^* for a E Next{K)) be the open and convex polygon with vertexes 
{xk,xl) {xk respectively) and {a(.,K,L)i<i<i- Observe that 



The discrete gradient V^m/i of a constant per control volume function Uh is defined as the 
constant per diamond Tk,l IR'-valued function with values 



Notice that : 

• The Z-dimensional mesure \Tk,l\ of Tk,l equals to j [ct^.l] dK,L- 

• The semi-norm coincides with the L'^{fl) norm of VhUh, in fact 

\Ur-a 



l|V.«.||i.(n) - E E / NnUul'dx^l'Y E 1^^.^ 



K\ 



KeTLeN(K) 
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• Let Fk^l for an arbitrary IR' vector associated with the interface aK,L satisfying F^^l = Fl,k- 

— * 

We denote by £h the set of interfaces aK,L- Then, a discrete field {FK,L)aK L&Sh assimilated 
to the piecewise constant vector function 

XI ^L,KXTk,l- 

The discrete divergence of the field F^ is defined as the discrete function = diVhFh with the 
entires 

diVKFh — 7^ 2^ WkmFk^l ■ Vk,l- 

' ' LeN{K) 

A key point of the analysis of the two-point finite volume schemes is the following kind of 
discrete duality property : 

Lemma 3.1. For all discrete function Wh on which is null on dfl, for all discrete field Fh 
on VL, 

^\K\ WKdivxPh = - X \Tk,l\ ^K,LWh ■ Fk,l- 
KeT TK,Le£h 



Proof. 



KgT KeTLeN{K) 



K\ WKdivxPh = XI X] Wk,l\ wkFk,l ■ rjK,L 

KeT LeN{K) 

KeT LeN{K) 

^~ol^ 1^ l\(^K,L\dK,Ll ^ VK,L-Fk,L 



2 ^ ^ I 

KeT LeN{K) 
- X \TK,L\^K,LWh-FK,L- 



l^KT 



□ 



Next, we approximate Mj(s)Vj9 ■ riK,L, {i = 1, 2.) by means of the values sk, sl and Pk,Pl 
that are available in the neighborhood of the interface (Tk,l- To do this, let us use some function 
Gi of {a,b,c) E R,^. The numerical convection flux functions Gi, Gi e C(1R^,1R), satisfies the 
following properties: 

(a) Gi{-, b, c) is non-decreasing for all 6, c e IR, 
and Gi{a, ■,€) is non- increasing for all a,c eJR; 

(b) Gi(a, a, c) = —Mi (a) c and G2{a, a, c) = M2(a) c for all a, c e IR; 

(c) Gi{a, b, c) = —Gi{b, a, — c) and there exists C > such that (3.1) 
\Gi{a, b,c)\ < C [\a\ + \b\) \c\ for all a,b,ceJR 

(d) There exists a constant mo such that 
{G2{a, b, c) — Gi{a, b, c))c > mo|cp for all a,b,c E JR. 
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Remark 3.1. Note that the assumptions (a), (b) and (c) are standard and they respectively 
ensure the maximum principle on saturation, the consistency of the numerical flux, and the 
conservation of the numerical flux on each interface. Moreover, the last assumption (d) will 
he used to obtain the estimate of discrete gradient of the pressure p. Practical examples of 
numerical convective flux functions can be found in JW^ . 

In our context, one possibility to construct the numerical flux Gi satisfying (13. ip is to split 
Mi in the non- decreasing part Mi^ and the non-increasing part Mj^; 



(3.2) 



M,^{z):= / (M/(s))+cis M,^{z):=- {M/{s))~ ds. 
Jo Jo 

Herein, = max(s, 0) and s~ = max(— s,0). Then we take 

G,(a, b; c) = c+ (M,t(a) + M,^(6)) - (m,^(6) + M,^(a)) . 

which leads to, 

'Gi{a,b; c) = -Ml (6) c+ - {-Mi{a)) q- 
G2{a,b; c) = M2{b) c+ - M^ia) c 

Note that the function s i— )■ Mi(s) is non- decreasing, and the function s M2{s) is non- 
increasing which lead to the monotony property of the function Gj. Furthermore, depending on 
the assumption (I^ on the total mobility we have, 

(G2{a, b, c) - Gi{a, 6, c))c = M{b)c+^ + M{a)c-'^ > moc^. (3.3) 

The next goal is to discretize the problem (I2.ip - (I2.2I) . We denote by I) an admissible 
discretization of Qt, which consists of an admissible mesh of Q, a time step At > 0, and a 
positive number N chosen as the smallest integer such that NAt > T. We set, 

r := nAt forn G {0, . . . , A^} 

dptl ■■= ^(^r' - Pl^') forn G {0, . . . , iV - 1} 



dK,L 



1 r^^' 



pTl ■■= .n+i .+1 / PiOdC fornG{0,...,iV-l} 



Pl -Pk 

fp'S ■■= l^j' I fp{x)dxdt fornG{0,...,iV-l} 



At\K\ 
1 

At\K\ 



K 

ftK'-=lvrn^\ I jji{^)dxdt fornG{0,...,iV-l} 



gx,L := / (g ■ VK,Ly d-f{x) = (g ■ riL,K) d-f{x) 
Jk/l Jk/l 

A finite volume scheme for the discretization of the problem (I2.ip - (l2.2p is giving by the 
following set of equations with unknowns P = {p^^)k&Ti ^ ^ [0; ^] ^-^d S = (s^^)i^g7-, n G 
[0, N], for all is: G r and n G [0, A^] 



Pk = J]^ J Po{x) dx, s?^ = ^ y So{x) dx, 



(3.4) 
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and 

L&N{K) 



E PtlGi{sl^\ dpp;l) + f/;/^) + \K\ p{pV-')st'rp'S = 0' (3-5) 

LeN{K) 



LeN(K) 



LeN{K) 

where the approximation of / p'^{p^~^^)Mi{s"'~^^)g ■ r]K,Ld'^{x) by an upwind scheme : 

JdK 

Kk= E Kk,l= E (p'(Pr')A^i(^r^)gi^,L-p^(p2^^)Mi(.2+^)g^,^), (3.7) 

L£N{K) L&N{K) 

and similarly F^^J^ the approximation of / P2^2(s"^^)g ■ r]K,Ld'^{x) such that 

4^^^'= E 4S= E (p2M2(.2-^l)gi.,L-p2M2(4+l)g^,^). (3.8) 

LeN(K) L&N{K) 

Note that the numerical fluxes to approach the gravity terms Fi , F2 are nondecreasing with 
respect to sk and nonincreasing with respect s^. 

We extend the mobility functions s ^ Mi{s) and s ^ M2{s) outside [0,1] by continues 
constant functions as follows, The approximate solutions, pst,h, sst,h '■ IR-^ x ^7 — )■ IR given for 
a\\K eT and n E [0, A^] by 

Pst,h{t, x) = pI+' and sst,h{t, x) = (3.9) 
for all X G -f^ and t G {nAt, (n + l)At). 
The main result of this paper is the following theorem. 

Theorem 3.1. Assume that (H^-(HE) hold. Let {po,so) G L'^{n,]R) x L°°{n,]R) and < 
So < 1 a.e. in Q. Then there exists an approximate solution {pst,h, sst,h) to the system liS. 5\) - 
l\3. 6\) . which converges (up to a subsequence) to {p,s) as {6t,h) — > (0,0), where {p,s) is a weak 
solution to the system /I2.1\) - l[2^) in the sense of the Definition \2.1[ 

4 A priori estimates 



We are now concerned with a uniform estimate on the discrete gradient of /3(s), and on the 
discrete gradient of the global pressure p. 
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4.1 Nonnegativity 

We aim to prove the following lemma which is a basis to the analysis that we are going to 
perform. 

Lemma 4.1. Let {s^x)KeT ^ [0, 1]- Then, the solution {s]^) Ki=T,n€{o,...,N} , of the finite volume 
scheme f l3.4p - fl3.6l) remains in [0,1]. 

Proof. Let us show by induction in n that for all G T, > 0. The claim is true for n = 
and for all K E T. We argue by induction that for all K E T, the claim is true up to order n. 
We consider the control volume K such that = min {s^'^^}^^^, and we seek that > 0. 
For the above mentioned purpose, multiply the equation in (13. 5p by —{s^^)^, we obtain 



LeN{K) ^'^ 
L£N(K) 

- FiT^'V-y - \K\ p{pT)sl^'fpAsl^')- = 0, (4.10) 
Observe that Pi^s"^^) — Pi^s"^^) > (recall that (3 is nondecreasing) . Which implies 



E V^(/3(^2+')-/3(^S:+'))(^x^')">0. (4.11) 

LSN{K) ^'^ 



The numerical flux Gi is nonincreasing with respect to s^"*"^ (see (a) in (13. ip ). and consistence 
(see (c) in (13. ip ). we get 

Using the identity s]^^ = {s]^^)~^ — {s]^^)~ , and the mobility Mi extended by zero on ] — oo, 0], 
then Mi(s^+i)(s^+i)~ = and 

- Fi:^'\sl^'r - \K\ p{pl^')sl^'f^pS{sl^')- 

= E p'(p2"'')^i(^2"'')gL,/r(4+')- + ii^iP(pr')/p,i^i(^r')-r>o. (4.13) 

L£N{K) 

Then, we deduce from (I4.10p that 

\K\ < 0, 



and from the nonnegativity of s^, we obtain (s^^) = 0. This implies that > and 

< s^+i < s2+i for all n G [0, - 1] and LeT. 
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To prove that < 1 for aWn G [0, — 1] andi^' G T. We argue by induction that for all 
For the mentioned claim, we multiply the equation in f l3.6[l by {s]^^ — 1) 



K E T, s]^ < 1. Let the control volume K such that = max {^^'''^j^^gT-, and let us show 



that s]+^ < 1 
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„n+l n I-. I 

L€N{K) 



+ |K| (4+1 - i)/p,i.(4^i - 1)+ = - jr^Ki^v-' - 1)^ (4.14) 



Since (3 is nondecreasing, we get /3{s2 ) — P{sk ) ^ 0- This implies 

- E ^(/3(^2^')-/3(4"''))(4"'^-i)^>o. (4.15) 

LgN{K) ' 

Next, we use the fact that the numerical flux G2 is nondecreasing with respect to and 
consistence (see (6) and (c) in (13. ip to deduce 

G,{sl^\ sT'- dpl^l) [st' - 1)+ > G,{sl^\ sl^'; dp^^l) {s^^' - 1)+ 

= dp^iM2(.ri)(.r'-l)+ = 0, ^ • ^ 

now, we rely on the extension of the mobility M2 by zero on [1, 00 [, thus M2{s'^'^) {s^^ — 1)"*" = 
0, to deduce 

4x'^(^r^-l)^= E P2M2{st')gKAsl^' - > (4.17) 

LeN(K) 

It is clear that the production source term in the left hand side of (]4.14p is nonnegative and 
the injection source term on the right hand side is nonpositive. 
Using the above estimates to deduce from (14.141) that. 



n+l _ n 



\K\ 



At 



((4+1 - l)(4+i - 1)+ - (4 - l)isT' - 1)+) < (4.18) 



Using again the identity (4+^ - 1) = (4+^ - 1)+ - (s^+i - 1)", and that 
< 1 to deduce from (I4.18P that (s^^^ — 1)+ = 0. Consequently, we obtain 

s2+i < s^+i < 1 for alln G [0, - 1] andL G T. 

□ 
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4.2 Discrete a priori estimates 

Let us recall the following two lemmas :: 
Lemma 4.2. (Discrete Poincare in equality JfJ 



Let n be an open bounded polygonal subset of M!^, d = 2 or 3, T an admissible finite volume 
mesh in the sense given in the section\^ and let u be a function which is constant on each cell 
K gT, that is, u{x) = uk if x E K, then 

where |H|//, (fj) is the discrete norm. 

Remark 4.1. (Dirichlet condition on part of the boundary) This lemma gives a discrete 
Poincare inequality for Dirichlet boundary conditions on the boundary dVL. In the case of 
Dirichlet condition on part of the boundary only, it is still possible to prove a discrete Poincare 
inequality provided that the polygonal bounded open set VL is connected. 

Lemma 4.3. (Discrete integration by parts formula) Let Vl be an open bounded polygonal subset 
ofMf'', T an admissible finite volume mesh in the sense given in the subsection\3i Let Fk/Li K G 
T and L G N{K) be a value in JR depends on K and L such that Fk/l = —Fl/k, and let ip be 
a function which is constant on each cell K eT, that is, ^{x) = ipx if x E K, then 

Fk/lVk = -\Y Fk/l{vl-Vk) (4.19) 

K&T L(^N{K) KeT L&N{K) 

Consequently, if Fk/l = CLK/iipL — bx), with Qk/l = o,l/k, then 

Y Y ^K/LipL-hK)VK = -]^Y Y ^K/LipL-hK){'^L- Vk) (4.20) 
KeT L&N{K) K&T LeN(K) 

Proof. The sum J^kgt '^l&n(k) reorganized by edge. In fact, on each edge ax^L between 

the mesh K and L, there are two contributions : from K to L named Fx^ifK and from L to 
K named Fl^k'^l, then 

Y Y Fkil^k = Y^^^I^'P^^^^I^^^^ ^^-21) 
Using now the fact that F^/l is antisymmetric, then we have 



Y Y ^K/L^K = Y ^K/l{^K - Vl) 

KeT LeN{K) (TK,L 



(4.22) 



2 



Finally, reorganise the last summation on edge by mesh, we obtain exaclty fl4.19p . The equality 
f !4.20p is a direct consequence of f l4.19p F^/l = dK/Li^L — bx) = o-K/hibK — &l) = ~Fl/k. □ 
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We derive in the next proposition, the main uniform estimates on the discrete gradient of 
the capillary term and the discrete gradient of the global pressure p. 

Proposition 4.1. Let {p^, s']^)KeT,ne{o,...,N}, be a solution of the finite volume scheme f l3.5p - 
(13. 6p . Then, there exist a constant C > 0, depending on Q, T, sq, po and a such that 



A) 



N-i I I (4.23) 

n=0 KeTLeN{K} ^'^ 



and 



KeT KeT 

1^-1 II (4-24) 

n=0 KeTLeN{K) 

where B'{s) = f3{s), and T-iij)) = g{p) + p{j>)p with g'{p) = —p{p). 

Proof. To prove the estimate (I4.23p . we multiply the gas discrete equation (13. 5p and the water 
discrete equation (13. 6p respectively by p^^^, g{p^^) = Tiip^^) — pip^^)p^^ and adding them, 
then summing the resulting equation over K and n, and this yields to. 



Ei + E2 + E3 + E4 + E5 = (4.25) 



where 



N-l 

^1 = E E 1^1 <Pk{{p{pI^')s-/' - p{pl)sl) pl^' + (.r' - si) g{pT')„ 

N-l I I 

E3 = E^^E E {pK':iGi{si^\sr-'-,dpi+i)pi+'+ 

n=0 K&T L£N{K) 

G,{s-jt\sl^\dp-+l) g{pl^')), 

N-l 

^4 = E E {f^kS + 9{pr)) ' 

n=0 KeT 

N-l 

^5 = E E (piPK^'^K^'fpi' pT + - 1)/pS 9{pt')^ 



n=0 K£TL£N{K) 



N-l 



n=0 K€T 



fit' 9{pT'' 
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To handle the first term of the equahty (14.251) . let us prove that : for all s > and > 0, 

{p{p)s - pip*)s*)p +{s- s^)('H(j9) - p{p)p) > nip)s - nip*)s\ (4.26) 
Indeed, denote g{p) = 1-L{p) — p{p)p then g'{p) = —p{p), 

{p{p)s - p{p*)s*)p +{s- s*)(n{p) - p{p)p) 

= s{H{p) - s\p{p*)p + g{p)) = sH{p) - s*H{p'') + s^(?/(p^) - p{p'')p - g{p)). 

We have to show that 

n{p*)-p{p*)p-9ip)>o. 

We expand this quantity as follows, 

n{p*) - pip*)p - g{p) = g{p*) + pip^'W -p)- g{p) = g{p*) - g{p) - g'{,p"W - P), 
as the function g is concave {g" (p) = —p'{p) < 0) , we get 

9{P) < 9{P*)+9'{P*){P-P*)- 
So, (14.261) is established, and this yields to 

E 1^1 'K-HiP^) - E 1^1 'K-Hip'K) < (4.27) 
Integrating by parts, see femma H73l we obtain 



N-l 

=0 KeTLeN(K) '^^'^ 



'pTM*' - Pin + iaipT') - 9(p'!t')))- 



Due to the correct choice of the density of the gas on each interface. 



_ i9{pl'') - 9{P1^')) 

/ 

we succeed to obtain, 



''''' ipV-'-Pl"') 



E2 = 0. (4.28) 

The choice of the density on the interfaces is the key point to vanish the dissipative term on 
saturation and obtain a uniform estimate on the discrete gradient of pressure p. 
Using the fact that the numerical fluxes Gi and G2 are conservative in the sense of (c) in (13. ip . 
we can apply lemma W^ a.iad we obtain 



Af-l 



n=0 KeT LeN{K) 
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Recall that inequality f l3.3p . 

'^2(0, b, c) - Gi(a, b, c))c = M(6)c+^ + M(a)c~^ > moc^ 
this with the hypothesis (H[5]) allow us to deduce that, 

N-l I I 

rnoP^E^^E E ^-^\pk^'-pI''\"<E3. (4.29) 

n=0 KeTLeN(K) 

To handle the other terms of the equality fl4.25p . firstly let us remark that the numerical flux 
satisfies F^~^^^ = —F^^f^^j^ and F2~^j^ = —F^^^, so we integrate by parts and we obtain 

n=0 K&T L&N{K) 

use now the fact that the mobilities and densities are bounded from (Hl2])-(Hl5]), and the map g 
is uniformly Lipschitz, we have, there exists a positive constant independent of At and h such 
that 

N-l 

n=0 KeT LeN{K) 

From the following inequality \o'k,l\ = (Ick.lMe'.l)^ 1°"^^!^ ^ g^j^^ apply the Cauchy-Schwarz 
inequality to obtain 

N-l 

\E^\ <Ci E E E WK,L\dK,L 
n=0 KeT LeN(K) 

N-l I I 

n=0 KeTLeN(K) 

N-l I I 

n=0 KeTLeN(K) ^'^ 

The last term will be absorbed by the dissipative term on global pressure from the estimate 

In order to estimate E5, using again the fact that the densities are bounded and the map g 
is sublinear {i.e.\g{p)\ < C\p\), we have 

\E,\ <c,Y,Atj2 \K\ {fp,K + w\ 

n=0 KeT 

we apply Holder inequality to deduce, 

N-l N-l 

\E,\<C,{Y,AtY,\K\\fp,K + f-%^\')-{Y,AtY,\K\ \pt^\'Y 

n=0 KeT n=0 KeT 

N-l 

n+l||2 \| 



<^^i(ii/p+Mi..(Q.)(E^^l|priL^(m) 



n=0 
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Now, from the discrete Poincare inegqality lemma leads to, 

7V-1 



\E,\<C,{J2^t\\pl^%J 



n=0 



where C2 is a constant depends only on \\fp + fi\\L2fn )■ Finally, under the assumption (Hl3]) 



(Qt)- 

on the source terms and as an application of Young's inequality {a ■ b < rjo? + |^.), we get 



N-l 



\E,\<C^+'^Y.^'Y: E ^^W-pT'\' (4-30) 

n=0 KeTLeN{K) 

this estimate fOOjl . with f OTj) . fOHj) and ( ICTj) achieve the proof of fOSj) . 

To prove the estimate (14.241) . we multiply the water discrete equation in (13. 6p by 
then summing the resulting equation over K and n, and this yields to 



where 



^1 + ^2 + ^3 = (4.31) 



= E E 1^1 <i''<('K^' - 'k) ■ Pisit'), 

n=0 KeT 

n=0 KeTLeN(K) 

N-1 

n=0 A'gT L(iN{K) 
N-1 

+ E E ■ /3(^r^^ 



K 

n=o i^er 



n=0 iCGT 

Let B{s) = / (3{r)dr. From the convexity of -B(s) (recall that /3"(s) = a(s) > 0), we 
Jo 



obtain 

N-l 



= E E 1^1 ('K^' - 'Kmst') 

n=0 K£T 

N-l 

> E E 1^1 (^('K^'^ - ^(^A-)) ^^-32) 
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|2 



Applying lemma lA^ we obtain 

n=0 KeTL€N{K) 

The other terms in the equahty fl4.3ip can be treated as fl4.30p . using Holder's and Young's 
inequalities with the help of lemma [^^2] and the assumptions on mobilities (Hl2]), source terms 
(Hl3]) and densities (HjS]) to get, 



N-l II 

\j,\<c+\Y.^tY: E ^l/3(4'^^)-/3(^nr (4.34) 

n=0 K£TL€N(K) ^'^ 



Now, collecting (021) . ( 103ll and (HIMD we obtain 



Af-1 



n=0 KeTLeN{K) 

for some constant C > 0. This concludes the proof of proposition 14. 1[ □ 



5 Existence of the finite volume scheme 

The existence of a solution to the finite volume scheme will be obtained with the help of the 
following lemma proved in [?] and [?]. 

Lemma 5.1. Let A be a finite dimensional Hilbert space with scalar product [■, ■] and norm 
\\-\\, and let V be a continuous mapping from A into itself such that 

[^(O,e]>0/or ||e||=r>0. 

Then there exists ^ E A with \\^\\ < r such that 

p(0 = 0. 

The existence for the finite volume scheme is given in 

Proposition 5.1. Let T> be an admissible discretization of Qt- Then the problem (I3.5p -( l3l6l) 

admits at least one solution {px, SK)(K,n)enRx{Q,...,N}- 
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Proof. At the beginning of the proof, we set the following notations; 

M := Card{T) 

We define the map 7^ : x IR^ — ^ IR^^ x IR-^, 

Th(sM,PM) = {{Ti,K}KeT, {T2,K}KeT) where, 

LeN{K) ^'^ 

+ E Pk:IG,{sI^\ si''; dpl^l) + F^lji) + \K\ p{pl^')sl^' f^^^', (5.35) 

LeN{K) 

„n+l _ n I I 

+ E G2{sl+\sl'';dpl^l) + F^:^^]; + \K\{sl+' (5.36) 

L£N{K) 

Note that Th is well defined as a continues function. Also we define the following homeomor- 
phism : X ^ X such that, 

J^iPM,SM) = ipM,VM) 

where, Vm = {gip]/') + Pisl+')}Ker- 

Now let us consider the following continues mapping Vh defined as 

VhiPM: Vm) ^ThO J^~'{pm: Vm) 
= Th{sM,PM)- 

Our goal now is to show that, 

[Ph{PM:VM): {PM:Vm)] > 0, for 1 1 (p^ , ^^^I ) 1 1 jj^2^ = T > 0, (5.37) 

and for a sufficiently large r. 
We observe that 



[Vh{pm,vm), {pm,vm)] E 1^1 'K^'Hipl^') - X7 E 1^1 'kH{pI) 

+7^El^l^(^^"^)-7^E^(^"-) 



At ^ ' ' ' ^ ' At 

KeT KeT 



+C{m,,pJ + 1/2 ||/3(.r)||^^(^) - C, 
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for some constants C^nio, pm),C > 0. This implies that 

[Ph{pm, vm), {PM, vm)] > - ^ E 1^1 'KHipl) - ^ E ^('k) 

Ker KeT (5.38) 

for some constants C,C' > 0. Finally using the fact that g is a Lipschitz function, then there 
exists a constant C > such that 

U{pl^'}KeTA9ipl^') + Pisl^')}KeT)\\j^^M < + • 

Using this to deduce from f l5.38p that fl5.37p holds for r large enough. Hence, we obtain the 
existence of at least one solution to the scheme (I3.5p -f l3?6|) . □ 



6 Space and time translation estimates 

In this section we derive estimates on differences of space and time translates of the function 
4'hPij>h)shB{sh) which imply that the sequence (j)hp{ph)shB{sh) is relatively compact in L^{Qt)- 
We replace the study of discrete functions U'^ = (l)hp{ph)shB{sh) (constant per cyhnder 
Qli '■— {i^,t^^^)] X K) by the study of functions tJ^ = (f)hp{,Ph)ShB{sh) piecewise continuous in 
t for all X, constant in x for all volume K, defined as 

U\t,x) = Y,Y. -^{{t - n^mi^' + ((n+l)At-t)t/]^) IlQ.(t,x). 

n=0 XeTh 

For a given discrete field J^h '■= J2aK l •^k,l^Tk discrete divergence is defined as a discrete 
function with entries on each control volume K] 



' ' L€N(K) 

Observe that we can write the discrete scheme fl3.5l) - fl3.6p in the following from: 

- CllV;,^^^^^ + 

— div/,y-2,/i + h,h ■ 



At 



where, 



n+l n+1 

o-K,L 

^ -(p^(pS:+^)M,(4+i)g;,,^-p2(p2+^)M,(.ri)g^,^)) ■r^^,aT,,, 
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^ ^aK,L 0-K,L 



0'K,L 

Ji,h ■— P\Pk J^K J P,K 

-pn+1 /ora+l 1 ^ \T<"\ 

J2,h ■— ^)JP,K \^\Jl,K- 

We also extend by the constant in time value U^'^h.+i [At{Nh + l) , +oo); as to J-'i^h, 
J^2,h, fi,h and f2,h, they are extended by zero values for t > At{Nh + l). The above definitions 
permit us to rewrite the equations fl3.5p - fl3.6l) under the form 

dt(j)hp{Ph)sh = diYhfi,h + fi,h ,^ 

(6.39) 

dt(t>hSh = di'^hJ^2,h + f2,h, 

where ^i^y ?2h-- fih /2/1 respectively the discrete functions of values vF"^^, ^2h^^ 

/{^+^ and /^f^ on each interval 

These equations are satisfied in iy^'^(lR''^) in time, for a.e. x & Q. 

Lemma 6.1. There exists positive a constant C > depending on Q, T, uq and Vq such that 



IL 



\U{t,x + y) - U{t,x)\'^ dxdt < C \y\ {\y\ + 2h), (6.40) 



'n'x{o,T) 

for all y G IR^ with Vl' = {x E fi, [x, x + y] C Vl}, and 



\U{t + T,x) - f/(t,x)|^ dxdt < C(r + At), (6.41) 



Qx(0,T-t) 



for allr G (0,T). 

Proof. The proof is similar to that found in, e.g, |8]. 

Proof of (I6.40p . First to simplify the notation, we write 

instead of 



Let ?/ G IR^ X G n\ and L G N{K). We set 



1, if the line segment [x, x + y] intersects aK,L, K and L, 
0, otherwise. 



y 

Next, the value Cg-^ , is defined by Co-^, r = -r~i ' Vk l with Co-,, , > 0. We observe that (see for 

\y\ 

more details [TU] ) 

(3^j,,^{x) dx < m{aK,L) \y\caKL^ 



(6.42) 

2^ l^<TK,L{^)(^'^K,LdK,L < \y\ + 2/i. 



Finite volume scheme for water gas 



21 



With this and an apphcation of the Cauchy-Schwarz inequahty leads to 

\U''{t,x + y)-U^{t,x)\^ dx 

{o,T)xn' 

<i\y\ + 2h) J2AtJ2 — 7^ / P^kA^) dx 

n=0 a^,, C^K,L<^K,L Jw (g43) 



< 



\y\i\y\+2h)J2^tJ2 



Ccrjf^dK,L 
n=0 CTK.L 



<C\y\i\y\+2h)J2^tY. 



Cajf]^dK,L 

for some constant C > 0. In addition, we have 

/ / \U^{t,x+dx)-U''{t,x)\dxdt <2 / / \U''{t,x+ dx)-U^{t,x)\ dxdt 
Jo Jn' Jo Jn' 

+2Ath [ \U^ix)\dx, 
Jn'^ 

where Uq = p{po)soB{so) and ^2'^ = {x E Q, dist(x,r2') < |A|}. By (16.401) . the assumption 
Ath — 7- as — 7- and the boundedness of (mo)/i in L^(fi^), then the space translates of on 
fl' are estimated uniformly for all sequence {hi)i convergent to zero. In the sequel, we drop the 
subscript i in the notation. 

Proof of fl6.4ip . Now we show a uniform estimate of the time translates of {U^)h '■ 

P + OO P 

for all A e (0,r], / \U''{t + A,x)-U''{t,x)\dxdt<u{T) (6.44) 
Jo Jn 

uniformly in h. Here u : — > M"*" is a modulus of continuity, i.e., lima;(r) = 0. 

Let us construct cD(-) verifying (I6.44p . First fix h and fix a G (0,r]. Denote by I'^{a) the 
left-hand side of (ICTD . For t > 0, set W''{t, ■) = U''{t + A, ■)-U''{t, ■). Notice that W''{t, ■) = 
for large t. 

Take a standard family (p<5)<5 of mollifiers on defined as ps{x) := 5^''p{x/5), where p is a 
Lipschitz continuous, nonnegative function supported in the unit ball of JR\ and J-^i p{x) dx = 
1. In particular, we have 

Here and throughout the proof, C will denote a generic constant independent of h and 6. For all 
t > 0, define the function (p{t, ■) : 1R' — > IR by (p{t) := ps * (sign t(;^(t)lLf7/). In order to lighten 
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the notation, we do not stress the dependence of on /i and 6. Discretize (p{t, ■) on the mesh 
Th by setting ipK(t) = ^ Jj^(p(t,x) dx; we denote (p^(t) the corresponding discrete function. 

Denote size (Th) '■= maxxgr^ diam (K). By the definition of ip(t,-), the discrete function ip^(t) 

is null on the set < x E 



dist (x, ^l') > 6 + size {%,) j, for all t. Thus for all sufficiently small 

h and 5, the support of f^(t) is included in some domain Q", Q" C Q. 
Note that 

Sit/'' = dt{p{ph)shB{sh)) = dt{p{ph)sh)B{sh) + p{Ph)shdtB{sh). 



Now for for all x G we multiply equation fl6.39p by \K\Lp{t)K-, integrate in t on [s, s + a], 
and make the summation over all K. Finally, we integrate the obtained equality in s over IR"*" 
to get 



J s 



Y,\K\^k{s)Wk{s) ds = 

J2 \K\ ipK{t)B{sK{t)) ( div,,[/l'(t)] + Uiit))K ) dtds 



K 

oo ps+A 



K 



(6.45) 



oo ps+A 



+ 1^ J Yl 1^1 VKit)piPkit))Skit)P{sKit)) i^dWKlf^it)] + if^it))Kj dtds. 

Denote by I^{a) the left-hand side of (K4^ . Denote Q" = (0, {Nh+l)At)xn" . Using hypothesis 
the definitions of discrete norms and the Fubini theorem, we get 



/'(a) < Ca 



7^ 



(max max I Vx l'p'' (t ) I + 



\VK,LB{sh{m) + ||¥^'1Uoo(Q")||riUHQ") 



Now the the L2^^([0,r] x fi) bounds on {P)hXf^)h, the bounds |</?(t, ■)! < 1 and |V(/?(t, ■)! < 
yield the estimate 

/^^(a) < Ca(1 + r'-^) (6.46) 
for all /i and 8 small enough, uniformly in h. Now, notice that by the definition of (pxit), 



\K\(^\WK{t)\-WK{t)y^K{t)) = \K\ \W\t,x)\ - WK{t)j^{t,x)dx 



K 



\W^{t,x)\-W^{t,x)^{t,x)]dx 



(6.47) 



therefore 



I\a) - I^{A) 



Jo. 



W^{t, x) I - W^{t, x)^{t, x) ) dxdt 



(6.48) 



Starting from this point, the argument of Kruzhkov Ref. [?] applies exactly as for the "con- 



tinuous" case. Set Uc := <x E JR 



dist (a;,9r]') < 6}: notice that Ui C Q" C Q for all 
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6 small enough. Notice that without loss of restriction, the boundary of Q' can be chosen 
regular enough so that to ensure that meas (f/^) goes to zero as 5 — )■ 0. By the result of 

Step 1 of the lemma and the Frechet-Kolmogorov theorem, the family (^J \W'^{t, ■)\dt^ is 

relatively compact in Ll^^{Q). Therefore these functions are equi-integrable on Q", so that 

P + OO P 

/ / \W^{t, x)\dxdt < uj{5) uniformly in h, with \ims^oUj{S) = 0. Now by the definition of 
Jo Ju'g 

(f, from formula ( I6.48P we deduce that 

\W\t,x)\dxdt 



\I'^{a)-I^{a)\<2 



+ 



JU: 



Jn'\u'g 



\W^{t,x)\-W''{t,x){ps* sign W'' (t) ) (x) 



dxdt. 



the first term in the right-hand side accounts for the action of the truncation Hqi in the definition 
of (fi. Using the standard properties of ps, we infer 



\I\a)-I^{a)\<2co{6) 



'0 Jn'\u', JWI 



\W''{t,x)\ - W''{t,x) signW''{t,y) 



dydxdt. 



Now note the key inequality: 

Va,6 G 

Setting a := {x — y)/6, we infer 



a — a sign 6 


< 2 1 







P + OO P P 

I\a)-I^{a)\ < 2u;{5) + 2 / / ps{x-y)\W\t,x)-W\t,y)\ dydxdt < 

Jo Jn' JlR 



+ 00 



(6.49) 



< 2u{6) + 2 / ^p{a) / / \u''{t,x)-u''{t,x-6a)\dxdtda < 2u{5) + 2u{5), 
JJR' Jo Jq' 



where u{-) is the modulus of continuity controlling the space translates of in Q'. Indeed, 
by Steps 1 and 2 of the proof, one can choose u{-) independent of h. Combining f l6.46p with 
( I6.49p . we conclude that the function 

u{t) := inf C { r (1 + S'^-^) + 2u}{6) + 2uj{6) } 



upper bounds the quantity I^. Because uj{t) tends to as r — )■ 0, this proves fl6.44p . □ 
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7 Convergence of the finite volume scheme 

Proposition 7.1. There exists a subsequences, still denoted {Uh,Sh)h, such that, as h 



Uh — > U strongly in L^{Qx) and a.e. in Qt for all p > 1, 
Sh — y s strongly in L'^{Qt) for all p > 1, 
"^hPish) V/3(s) weakly m {L\Qt))\ 
^hPh — > weakly m {L^{QT)f, 
Ph — y p weakly in L?{Qt)- 



Furthermore, 



Finally, we have, 

fl{Ph)f2{Sh) 



Sh — > s a.e. in Qt, and < s < 1 a.e. in Qt, 
U = (f)p{p)sB{s) a.e. in Qt 



f\{p)f2{s) a.e. in Qr.V/iJa e C°(]R) such that /2(0) = 0. 



(7.50) 
(7.51) 
(7.52) 
(7.53) 
(7.54) 
(7.55) 
(7.56) 



(7.57) 
(7.58) 



(7.59) 



Proof. For the first convergence f l7.50p it is useful to introduce the following inequality, for all 
a,6 e IR, 



/ \ea + {i-e)h\de>\{\a\ + \h\) 

Jo ^ 



Applying this inequahty to a = u'^^^'' — u^, b 



/ / \U''{t,x) -U^{t,x)\dxdt<2 / / \U^{t+Ath,x) -U^{t,x)\dxdt. 

Jo Jw Jo Jq' 

Since Ath tends to zero as /i — )■ 0, estimate ( I6.44p in Lemma 16.11 implies that the right-hand 
side of the above inequality converges to zero as h tends to zero, and this established fl7.50p . 

By the Riesz-Frechet-Kolmogorov compactness criterion, the relative compactness of {U^)h 
in L}{Qt) is a consequence of Lemma [6.11 Now, the convergence fl7.5ip in L^{Qt) and a.e in 
Qt becomes a consequence of fl7.50p . Due to the fact that is bounded, we establish the 
convergence in L}{Qt)- 

In order to prove the third convergence f l7.52p . we reproduce the previous lemma EH] for 
= ShB{sh), and as an application of the Riesz-Frechet-Kolmogorov compactness criterion 
we establish f l7.52p . 

For the weak convergence of the discrete gradient of the global pressure, let us recall the 
piecewise approximation VhPh of Vph in Qt- 



(n-l) 



from the definition of we deduce 



VhPh{t,x) 



pn pn _, 

I^^^^Vk,l if {t, x) G (r, r+1) X Ti 



K,L, 

if {t,x) G (r,r+i) X Ti?;;, 
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for a\\ K E T and < n < Nh- It follows from proposition 14.11 that, the sequence (VhPhjh is 
bounded in {L^{Qt))^, and as a consequence of the discrete Poincar inequality, the sequence 
{ph)h is bounded in L'^{Qt)- Therefore there exist two functions p G L'^{Qt) and ^ G (L^(Qt))^ 
such that f l7.55p holds and 

^hPh — > C weakly in {L^{QT)f. 

It remains to identify Vp by C, in the sense of distributions. For that, it is enough to show as 
/i ^ 0: 



Eh-= / VhPh-^dxdt+ / / phdivipdxdt — ^ 0, V(/7 G ^'(gT) • 

Let h be small enough such that ip vanishes in T^); for all K eT. In view of riK,L = K 
we obtain for all t G (t", 

/ divy9(t, x) = / div x) 

= XI 5Z / ^it,s) ■7]K,Lds 

= E (Pk-pD [ ^{t,s)-VK,Lds. 

T ^ AT/ J CTK J, 



K&T LeN(K) 

Now, from the definition of the discrete gradient 



Then, 



I VhPh^{t,x)dx = ]^^ X / VhPhf{t,x)dx 

= IY1 Y1 -T—(pI-Pk)[ ^{t, ■ Vk,l dx 

^^^=0X1 ^kApI-p1)( / vit,s)-r]K,Lds 



I V{t,x) ■ VK,Ldx 

^K,lJTk,l 

Due to the smoothness of (p, one gets 

^ ' <C h, 



[ <^{t,s) ■T]K,Lds - ^ [ ip{t,x) ■ T]K,Ldx 

and the Cauchy-Scharwz inequality with proposition 14.11 

\Eu\<chY,^tY. E M\pi-pi\ 

n=0 K&T L£N{K) 

<ChY,^tYl E \''K,L\dK,L 
n=0 KeT LeN{K) 

< Ch\Q\T. 
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Now, for the identification of the hmit fl7.58p : 
Due to the monotonicity of the function p, we have 



kp(,Ph)shB{sh) - (f>hP{v)shB{sh)) dxdt > 0, Wv e L^{Qt) 



this with the strong convergence (17.511) and the weak convergence (17.551) lead to, 



{U - <pp{v)sB{s))dxdt>Q, WveL^Qt) 



Q 



T 



Finally, choose v = p + 6w with 6 g]0, 1] and w G L'^{Qt), then 



(U - 0p(p + 5w)sB{s))w dxdt > 



letting 6 goes to zero, we establish (17.581) . 

To conclude the a.e. convergence (17.591) .on one hand, when s = a.e., fi{ph)f2{sh) 

= /i(p)/2(s) a.e. (since /2(0) = and /i is bounded). On the other hand, when s/i — > s 7^ 0, 
in light of (17.511) we have fijp h)^ fijpi) a.e.. Then, fi{ph)f2{sh) /i(p)/2(s) since /i, /2 are 
continuous and this establish ( I7.59I) .D 



Theorem 7.1. Assume (Hl\)-(I{5\) hold. Then the functions p,s defined in proposition 7.1 
constitute a weak solution of the system l^2.1\) - [27^) . 



Proof. Let T be a fixed positive constant and ip G I){[0,T) x Q). 

• Convergence of the discrete water equation 
For the discrete water equation, we multiply the equation ( 13. 6p by Atip(t'^~^^, xk) for all G T 
and n G {0, . . . , A^}. This yields 



(T^ + + (T^ + (T^ + (T^ + (T^ = 
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where 



N-l 

^1 = E E 1^1 M^K^' - siMr+\xK), 

n=0 K£T 



N-l 

n=0 KgTL€N{K) ^'-^ 

N-l 

n=0 KeTLeN{K) 

N-l 

n=0 KeT LeN{K) ''^1^ 

LeN{K) -^^1^ 

N-l 



n=0 K^T 
N-l 



n=0 KeT 

Performing integration by parts and keeping in mind that (fi{T, xk) — for all X e T, we 
obtain 

Af-l 

= - J] J] cPKS''+\ip{r+\xK) - V{t\XK)) - Y 1^1 <Pks'MO,Xk) 

n=0 K&T K&T 

N-l „tn+i „ „ 

= / / (l)KSKdt'p(t,XK) dxdt - Y / 0ii:So(a;)</'(O,xx) rfa:; 

— ■ ^1,1 *~'l,2- 

Let us also introduce 



•^i'* ^ ~ E E / / <t>Ksl:dtip{t,x)dxdt - I (f)hSo{x)(p{0,x) dx 
Now and due to the fact that the saturation and the porosity functions are bounded, we have 

N-l rt"+^ 



Y E ^^^K / / (dt(p{t,XK) - dt(p{t,x)j dx 

N-l tn+l 

<0iEE / / dt(p{t,XK) - dt(p{t,x 
n=o xer-^*" 
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using that the function (f is regular enough, we get 



hm 



^1,1 ^1,1 



Similarly 



xk) — V^(0, x)) dx 



(phSo{x){Lp{0, xk) - V5(0, x)) dx. 

By the regularity of there exists a positive constant C such that 
1^9(0, Xk) — V'(0, a;) I <Ch. This implies 



^1,2-^1$ <Ch(PiY^ I so{x)dx. 



Sending /i — )■ in the above inequality, we get 

^1,2 ^1,2 

Combining (17.601) with f l7.6ip . we obtain 



lim 



0. 



lim 

h-^0 



but, Cj^'* can be written equivalently. 



(f)hShdt(p{t, x) dx dt — / (j)hS^(p{0, x) dx 



(7.60) 



(7.61) 



(7.62) 



Since the bounded functions (j)h and converge almost everywhere respectively to (f) and s, 
and as a consequence of Lebesgue dominated convergence theorem, we get 



lim a = lim 



(f)sdt(p{t, x) dx dt — / (j)s'^(p{0, x) dx 



Now, let us show that 



Integrating by parts 



Af-l 



lim d 



Jn 



V/3(s) ■ Vip dxdt. 



(7.63) 



/3(s2+^) - /3(s^+^) v{t-+\xL) - ^{r^\xK) 



n=0 A'er LeN{K) 
N-l 



d 



K,L 



d 



K,L 



^ E E E l^^.^l ■ VK,lV^{r^\ Xk,l) ■ VK,L, 



n=0 KeT LeN{K) 
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where xk,l = (^xk + (1 — (^)xl, < 6* < 1, is some point on the segment ]xk,xl[- Recall that 
the value of Vk,l is directed by rj^^L, so 

since each term corresponding to the diamond Tk,l appears twice, 



6^2 = / / ^hl3{sh) ■ {V^)h dxdt, 
Jo Jn 



'0 Jn 
where 

Observe that from the continuity of ip we get (V(/9)/i V(/9 in L°°{Qt)- Hence the convergence 

f l7.63p is a consequence of fl7.53p . 

Now, we show the convergence of the flux, 

lim€.^ = - [ [ M2(s) Vj9 ■ Vip dxdt. (7.64) 
Jo Jn 

Perform integration by parts (14.191) . thanks to the consistency of the fluxes, we obtain 

N-l 

€l = --J2^tJ2 E G,isl^\sl'';dpl^l)[^{e-'\x,)-^{e^\xK)). 

n=0 KeT LeN{K) 

For each couple of neighbours K, L we denote s]^l the minimum of and s^^^ and we 
introduce, 

it'/ = --J2^tYl E M,{si^i)dpi^i{y,{r^\x,)-v{r^\xK)). 

n=0 KeT L£N{K) 

Define Sh and by 

Sh\{t",t"+^xTK,L ■- ^^i^K ■:^L h Sy^\{t",t^+^xTK,L ■- ^™-\^K ■:^L S' 

Now, £3'* can be written under the following continues form. 



^3'* = - / / M2{s^)VhPh ■ (Vv^)ft dxdt. 
Jo Jn 



'0 Jn 

By the monotonicity of /3 and thanks to the estimate (I4.24p . we have 

I I \p{-SH)-P{s^)?dxdt<Y^AtY^ E \tk,l\{(^k+') - (^{^i^': 



n=0 KgT L&N{K) 

N-l I I 



n=0 KeTLeN{K) ^'^ 



< Ch 



2 
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Since (3~^ is continues, we deduce up to a subsequence, 

I'^h a.e. on Qt- (7.65) 

Moreover, we have; s^^ < < 'Sh and Sh s a.e. on Qt- Consequently, and due to the continu- 
ity of the mobihty function M2 we have M(sfJ — > M(s) a.e. on Qt and in L'p{Qt) for p < +00. 
Using proposition 17. II f l7.54p and the strong convergence of (Vv?)/i to Vy?, we obtain that 

lim = - / / M2{s)Vp ■ Vip dxdt. 

It remains to show that 



Jn 



\im\(t^-€^'*\=0 (7.66) 

ft— >o 



By the properties of the numerical flux function (13.1 p we have 



Consequently, 



Jo Jn 



Applying the Cauchy-Schwarz inequality, and thanks to the uniform bound on VhPh and the 
convergence (I7.65p . we establish (I7.66p . Now, we treat the convergence of the gravity term 

limC^ = -/" [ p2M2is)g-Vifidxdt. (7.67) 
Jo Jn 

Perform integration by parts (I4.19p . 

Af-l 

n=0 K&T LeN{K) 

N-1 



XK. 



n=0 KeT LeN{K) 



We introduce, 



1 r 

(^ipir^\xL)~^ie^\xK)) 



Finite volume scheme for water gas 



31 



where s'^l := mm{s^+\ sV'^}- We have 

N-l „ 
n=0 K&TL£N{K) "^^'/^ 

((^(r+\xi)-</.(r+\x^)) 

relying on the assumption (Hl2]), that the mobihty functions are Lipschitz, and is a Holder 
function, we deduce that 

Jk/l 

= {m2{sT') - M,{sl^l)) [ (g • r^K^L^ d^{x) 
^ ^ Jk/l 

- (m,{sI+') - M,{sl^l)) [ (g . vk,l)- d^ix) 

^ ^ JKIL 

< C\g\\(JKr 



n+l _ n+1 



< C|g||o-i^,L| P{S 



this yields to 



N-l 



\€!l-€!l^*\<cY,^tY. H WK,L\dK,L\fi{sl^')-fi{st')\\VHiPh\. 

n=0 KeT LeN{K) 

Using estimate fl4.24p . and Cauchy Schwarz inequality, then 

1^4 - ^4*\ when 
For the convergence of the source terms, £5 + £g can be written equivalently. 



N-l 



^5 + ^6 = E E / / - ^)fpit,xMt-^\xK) 

N-l 

Now, we introduce 



dxdt 



n=0 K(^T 



fi{t, x)(p{t"'~^ \xk) dxdt 



N-l 



n=0 KeT 



-.n+l 
'K 



l)fp(t, x)ip(t, x) dxdt 



N-l 



n=0 K€T 



t" J K 



fi(t, x)ip(t, x) dxdt 
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Due to regularity of ip we obtain, 



^5 ' ^6 ^5 



<C(At + /.)(||/p|Ul(Q,)+||MUl(Q,)) 



and this ensures, 



We can write equivalently. 



lim 



^5 + 



{sh — l)fp{t, x)(p{t, x) dxdt + / fj{t, x)ip{t, x) dxdt 



Finally, by the convergence of the saturation function we get. 



{s — l)fp{t, x)(p{t, x) dxdt + / fi{t, x)(p{t, x) dxdt 

Qt JQt 



Convergence of the discrete gas equation 



Multiplying the discrete gas equation fl3.5p by At(y9(t"+^, a;^) for all K ^ T and n G 
{0, . . . , N}. Summing the result over K and n yields 



where 



+ + + Cl + = 0, 



N-l 



n=0 KeT 

N-1 



=0 KeTLeN{K) 



N-l 



C3 = E^^E E pS2Gi(4+\.ri;rfpg-2)^(r+\x;,j, 

n=0 KeT LeN{K) 

N-l „ 



n=0 iCGT 



L&N{K) 



- p\pT')m,{sI^') J2 [ {g-VK,L)-d^{x))^{r+\xK), 



L€N{K) 



n=0 KeT 
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Performing integration by parts and keeping in mind that (p{T, Xk) = for all i^T G T, we 
obtain 



N-1 



EE / / <pKp{pV-')s'}t'dMt,XK)dxdt 

E / (t>Kp{pV)A^{^^^K)dx 



nh nh 



Let US also introduce 



EE / / MPlt')sl^'dMt.^)dxdt 

E / 'pKp{pV)S^K^{Q^x)dx 



ph,* ph.* 
'-'1,1 ^ '-'1. 2 ■ 



Then 



Cl,2-Cf;2 = E / <PKP{PI)SI{^{0,XK)-V{0,x))dx. 



By the regularity of cp, there exists a positive constant C such that 
1^9(0, Xi^) — ip{0, x)\ < C h. This implies 



ph ph,* 

'-'1,2 ~ '-'1,2 

Sending — )• in the above inequality, we get 



< Ch\Q\ 



lim 



ph ph,* 

'-'1,2 '-'1,2 



0. 



(7.68) 



Now, due to the fact that the saturation function is bounded and the assumptions (I|T]),(Hl5])on 
the porosity and the density, we have 



ph ph,* 

'-'1,1 '-'1,1 



N-l „ 

E E MPt')si^' / 



n=0 K£T 
N-l 



K 



< 



pM E E 

n=0 K&T 



n + l 



K 



d(p{t, Xk) — d(p{t, x) ) dx dt 
dip{t,XK) — d(p{t,x) dxdt 
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using that the function ip is regular enough, we get 

hm Cl - Cf;; = 0. (7.69) 
Combining fl7.69p with f l7.68p . we obtain 

C^^-Cf'* =0 (7.70) 



hm 



but, Cf'* can be written equivalently. 



Jqt Jn 



Ci* = / (j)hpiPh)shdt(p{t,x)dxdt - / (j)hp{Ph)Sf,(p{0,x)dx. 



Since (l)hP{Ph)sh and (phP{Ph)^h converge almost everywhere respectively to (l)p{p)s and (l)pijP)s^, 
and as a consequence of Lebesgue dominated convergence theorem, we get 



h 

Now, let us show that 



lim = lim •* = / 

h,^o h^G Ja, 



(j)p{p)sdt(p{t, x) dx dt — I (j)p{p'^)s'^(p{0, x) dx. 



limC2^= [ [ p{p)V(5{s)-V^dxdt. (7.71) 
Jo Jn 

Integrating by parts 



N-l 



=0 KeTLeNiK) 

N-l 

2 



^-1 

oE^^E E \TK,L\pV:lVK,LPK''')-VK,N^it--'\xK,L)-VK,L, 



=0 i^GT LeN{K) 

where xk,l = Oxk + (1 — 6')a;L, < 6* < 1, is some point on the segment ]xktXl[- Recall that 
the value of V k,l is directed by ?7i^,L, so 

since each term corresponding to the diamond Tk,l appears twice, 

C2= [ I p{Ph)VhP{sh) ■ {V^)h dxdt, (7.72) 
Jo Jn 

where 
Define 

N-l I I 

n=0 KeTLeN{K) ^'^ 
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where = (piP^^) + P(^r^^)/2. We have : 

= (/3(^r^)p(^r^) - /3(4^^)p(i'r^)) 
+ Wist') + /3(4^^))(p(^r^) - p(^r'))/2. 

Then, can be rewritten 

= dI + 

where 

Af-l 

^3 = E^^ E ^^.^(/3(^r')p(pr')-/3(^r')p(px')) 

n=0 o-K,L&E 

ri=0 a-K,L&E 

where /J^"^/ = (/^(s^'^^) + /3(s^"^^))/2, recall that t^.l = Follow the same lines as in 

dlZlD, 

^3 = / /" ^h{p{sh)(3{sh)) ■ iV^)h dxdt, 
Jo Jn 

Di= [ [ Msh)Vhp{Ph) ■ {Vv)h dxdt 
Jo Jn 

Using dOSD and (gJH), we have 

Vh(p(ph)/3(s/,)) ^ V(p(s)/3(s)) weakly iriL^Qr), 
and using the fact that (Vip)h converges strongly in L'^{Qt), we have 

— ^ I I V(p(s)/3(s)) ■ dxdt. 
Jo Jn 

In order to handle the convergence of we are going to show 

Msh) ^ Pis) strongly in inL^Qr), (7.74) 

and 

^hP{Ph) — > Vp* weakly in inL'^^Qr)- (7.75) 
The sequence {p{ph))h is bounded then 

p{Ph) P* weakly in L^Qt). (7.76) 
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Using the fact that p'{.) is bounded, we have 

N-l 

N-l 

and using the estimate ( 17.54p . we deduce that Vhp{ph) is bounded in LF'^Qt) and converges 
weakly to a function ^ in -L^(fi); and from (17.761) we deduce that ^ = Vp* weakly. 
Recall that 

N-l 

f^i^h) = E I]^(^^"'')l^x]*"'*"+^1 '^(^) Strongly in L\Qt). 

n=0 K 

Let us show for ^(S/,) = J2n=0 lSE l3K^L^TK,L><]t",t-+^, 

Where = that 

Msh) - (3{sh) strongly in L\Qt) 

In fact. 



Sh) - ^{Sh)\\L^Qr) 
N-l 

= E^^ E (iTA^Lnif|i/3]^+2-/3rT+iTx,LnL|i/3]5,^-/3rip) 

n=0 CTK^L^E 
N-l 

n=0 (TK,L&E 

and from estimate fl4.24p . there exists a positive constant C such that 
which establish the desired limit. Then, the convergences fl7.75p . leads to 



limD^= [ [ (3{s)Vp* -Vip dxdt. (7.78) 
Jo Jn 



'0 Jn 



Finally, let us show C2 — D2 ^ 0. We have 

N-l 



n=0 aK,L 
~n+l ^"+1 tt'^+I 



where = p^^ — p^^ . This expression can be rewritten as 

C2 - ^2 = / ~p{Ph)Vhfi{sh) ■ (V0), dxdt (7.79) 
Jqt 
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Let us show that p{ph) — ^ strongly in L^{Qt)- 
We have 



and from hypothesis (HS]), the function p is monotone and uniformly Lipschitz, then there 
exists a positive such that 

pV:1<c\pi''-p"^'\- 

So, 

llp(p^)lli^(Q.) = E^^ E I^^.^IIpXlT (7-80) 

71=0 crK,L&E 
N-1 

<Cj2At ITkMV-'-pI^'I'), (7.81) 

using (14.231) . we deduce 

\\P{PH)\\l.^Q,)<Ch\ 

which goes to zero when h goes to zero. This convergence combined with the weak convergence 
fl7.54p and the strong convergence of (V0)/i in L°°{Qt) shows that 

Cl^ -D^ — ^ 0, when h^Q. 

□ 

Remark 7.1. In the case where K(x), the permeability tensor of the porous medium at a point 
X, considered to be 

K(x) = k{x)Xd 

where k is a scalar bounded function of the space, k{x) > k^ > and is the identity matrix. 
The main part is the approximation of the dissipative terms (capillary terms) on each interface 
o'K^L as follows: 
Denote by, 

1 f 

kx = T-fTT / k{x) dx. 
\K\ 

Now, we consider the following approximation 



K 



k{x)p{p)VP{s) ■ 7]K,Ldl ^ d*KLPK,J^{[5{sL) - P{sk)) 

cIk,l 



/ k{x)p{p)M,{s)Vp -mdl^ d*K^LPK,L{ - Mi{sL){dpK,LV + Mi{sK){dpK,Ly) 
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J OK J, 



where, 



and, 



d{XK,CrK,L)kK,L + d{XL,aK,L)kL,K 

dpK,L = ^-^r^(PL-PK) = {dpK,LV - {dpK,L)~, 
CL-- - 



K,L 



Sk,l-= / {g-'nK,L) da= / (g-?7L,ic) d'j{x) 

IK/L JK/L 

1 



Pk,l = \Pl~Pk Jpji 

.p{Pk) otherwise 

In fl7.82p ■ we take the harmonic average on the interfaces in order to ensure the conservation 
of numerical fluxes. 
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